Notes on adding a leapsecond to the leapsec_mmddyy.fits table

This notebook explains how to add a leap second to the leapsecond.fits file on the rare occasions that a leap second is added

References:

http://tycho.usno.navy.mil/leapsec.html
http://maia.usno.navy.mil/ser7/tai-utc.dat
https://hpiers.obspm.fr/iers/bul/bulc/
http://heasarcdev.gsfc.nasa.gov/howto/Leap_Second_Updates.html - Ed's instructions on who does what

See the pycaldb routine update_caldb_leapsec.py for the implementation of the code below



In [1]:

    
from astropy.io import fits as pyfits
from astropy.time import Time
from astropy.table import Table
import ftputil
from datetime import datetime
import time

Typical announcement of a Leap Second from https://hpiers.obspm.fr/iers/bul/bulc/: SERVICE INTERNATIONAL DE LA ROTATION TERRESTRE ET DES SYSTEMES DE REFERENCE SERVICE DE LA ROTATION TERRESTRE DE L'IERS OBSERVATOIRE DE PARIS 61, Av. de l'Observatoire 75014 PARIS (France) Tel. : +33 1 40 51 23 35 e-mail : services.iers@obspm.fr http://hpiers.obspm.fr/eop-pc Paris, 6 July 2016 Bulletin C 52 To authorities responsible for the measurement and distribution of time UTC TIME STEP on the 1st of January 2017 A positive leap second will be introduced at the end of December 2016. The sequence of dates of the UTC second markers will be: 2016 December 31, 23h 59m 59s 2016 December 31, 23h 59m 60s 2017 January 1, 0h 0m 0s The difference between UTC and the International Atomic Time TAI is: from 2015 July 1, 0h UTC, to 2017 January 1 0h UTC : UTC-TAI = - 36s from 2017 January 1, 0h UTC, until further notice : UTC-TAI = - 37s Leap seconds can be introduced in UTC at the end of the months of December or June, depending on the evolution of UT1-TAI. Bulletin C is mailed every six months, either to announce a time step in UTC or to confirm that there will be no time step at the next possible date. Christian Bizouard Head Earth Orientation Center of IERS Observatoire de Paris, France



In [2]:

    
NewLeapSecondDate = '2017-01-01T00:00:00' # in isot format
NewLeapSecond = 1.0



In [3]:

    
#
# this block retrieves the lastest leapsecond file from /FTP/caldb/data/gen/bcf
# based on the leapsecond file naming convention, "leapsec_mmddyy.fits"
#
LSdir="FTP/caldb/data/gen/bcf/"
host=ftputil.FTPHost('heasarc.gsfc.nasa.gov',"anonymous","mcorcoran@usra.edu")
genbcf=host.listdir(LSdir) # get directory listing
LeapsecFileList= [f for f in genbcf if 'leapsec' in f]
LeapsecFileYear = [y.split("_")[1].split(".fits")[0][4:6] for y in LeapsecFileList]
LeapsecFileYear = [('19'+y if int(y) > 50 else '20'+y) for y in LeapsecFileYear]
LeapsecFileMonth = [m.split("_")[1].split(".fits")[0][2:4] for m in LeapsecFileList]
LeapsecFileDay = [d.split("_")[1].split(".fits")[0][0:2] for d in LeapsecFileList]
maxjd=0.0
for i in arange(len(LeapsecFileList)):
    tiso=LeapsecFileYear[i]+"-"+LeapsecFileMonth[i]+"-"+LeapsecFileDay[i]
    fjd = Time(tiso).jd
    if fjd > maxjd:
        maxjd = fjd
        LatestLSF = LeapsecFileList[i]
print("Latest Leapsecond File = {0}".format(LatestLSF))









    



Latest Leapsecond File = leapsec_010715.fits



In [4]:

    
hdu=pyfits.open("http://heasarc.gsfc.nasa.gov/"+LSdir+"/"+LatestLSF)



In [5]:

    
tbdata=hdu[1].data



In [7]:

    
lsdate = tbdata['DATE']
lstime = tbdata['TIME']
lsmjd = tbdata['MJD'] # this is the MJD corresponding to the tabulated DATE and TIME
lssec = tbdata['SECONDS'] # elapsed seconds from MJDREF to MJD, with 1 day  = 86400 seconds
lsleapsec = tbdata['LEAPSECS']

# lsmjd[0] - Time((lsdate[0]+'T'+lstime[0]).strip(), format='isot').mjd



In [8]:

    
test=np.asarray(lsdate)
test=np.append(lsdate,'2017')
test









    Out[8]:





array(['1972-01-01', '1972-07-01', '1973-01-01', '1974-01-01',
       '1975-01-01', '1976-01-01', '1977-01-01', '1978-01-01',
       '1979-01-01', '1980-01-01', '1981-07-01', '1982-07-01',
       '1983-07-01', '1985-07-01', '1988-01-01', '1990-01-01',
       '1991-01-01', '1992-07-01', '1993-07-01', '1994-07-01',
       '1996-01-01', '1997-07-01', '1999-01-01', '2006-01-01',
       '2009-01-01', '2012-07-01', '2015-07-01', '2017'], 
      dtype='|S10')



In [9]:

    
lssec









    Out[9]:





array([  6.30720000e+07,   7.87968010e+07,   9.46944020e+07,
         1.26230403e+08,   1.57766404e+08,   1.89302405e+08,
         2.20924806e+08,   2.52460807e+08,   2.83996808e+08,
         3.15532809e+08,   3.62793610e+08,   3.94329611e+08,
         4.25865612e+08,   4.89024013e+08,   5.67993614e+08,
         6.31152015e+08,   6.62688016e+08,   7.09948817e+08,
         7.41484818e+08,   7.73020819e+08,   8.20454420e+08,
         8.67715221e+08,   9.15148822e+08,   1.13607362e+09,
         1.23076802e+09,   1.34110082e+09,   1.43570880e+09])



In [10]:

    
mjdref = hdu[1].header['MJDREF']
sec = (lsmjd[26] - mjdref)*86400 
print("{0:10.8e}".format(sec))









    



1.43570880e+09



In [11]:

    
outfile='leapsec_'+(time.strftime("%d%m")+(time.strftime("%Y"))[2:])+'.fits'
outfile









    Out[11]:





'leapsec_130716.fits'



In [12]:

    
time.strftime('%Y-%m-%d %H:%M:%S')









    Out[12]:





'2016-07-13 16:27:35'



In [13]:

    
tbdata.columns









    Out[13]:





ColDefs(
    name = 'DATE'; format = '10A'
    name = 'TIME'; format = '15A'
    name = 'MJD'; format = 'D'
    name = 'SECONDS'; format = 'D'; unit = 's'
    name = 'LEAPSECS'; format = 'D'; unit = 's'
)



In [14]:

    
#
# create new row for new leapsecond information
#
newdate = NewLeapSecondDate.split('T')[0]
newtime = NewLeapSecondDate.split('T')[1]
newmjd = Time(NewLeapSecondDate, format='isot').mjd
newsecs = (newmjd - mjdref)*86400
newLS = NewLeapSecond
NewLeapSecondRow=[newdate,newtime,newmjd,newsecs,newLS]
NewLeapSecondRow









    Out[14]:





['2017-01-01', '00:00:00', 57754.0, 1483228800.0, 1.0]



In [15]:

    
#
# append new leapsecond info to table and write to output file
#
tbdata=hdu[1].data
t=Table(tbdata) # convert hdu data to a python Table to add new row
t.add_row(NewLeapSecondRow) # add row of data
hdunew=pyfits.table_to_hdu(t) # convert table back to hdu (with minimal header)
hdunew.header=hdu[1].header # use header from original file
hdunew.header['COMMENT']='added 2016-12-31 leap sec by MFC'
hdunew.header['HISTORY']='updated by MFC 2016-07-12'
pyfits.writeto(outfile,hdunew.data,hdunew.header, clobber=True)



In [ ]:

    
#
# append new leapsecond info to table using FITS columns and write to output file
#
tbdata=hdu[1].data
t=Table(tbdata) # convert hdu data to a python Table to add new row
newdate = np.append(lsdate, NewLeapSecondDate.split('T')[0])
newtime = np.append(lstime, NewLeapSecondDate.split('T')[1])
newmjd = np.append(lsmjd,Time(NewLeapSecondDate, format='isot').mjd)
newsecs = np.append(lssec,newmjd - mjdref)*86400)
newLS = np.append(lsleapsec, NewLeapSecond)

t.add_row(NewLeapSecondRow) # add row of data
hdunew=pyfits.table_to_hdu(t) # convert table back to hdu (with minimal header)
hdunew.header=hdu[1].header # use header from original file
hdunew.header['COMMENT']='added 2016-12-31 leap sec by MFC'
hdunew.header['HISTORY']='updated by MFC 2016-07-12'
pyfits.writeto(outfile,hdunew.data,hdunew.header, clobber=True)



In [16]:

    
tbdata['DATE']









    Out[16]:





chararray(['1972-01-01', '1972-07-01', '1973-01-01', '1974-01-01',
       '1975-01-01', '1976-01-01', '1977-01-01', '1978-01-01',
       '1979-01-01', '1980-01-01', '1981-07-01', '1982-07-01',
       '1983-07-01', '1985-07-01', '1988-01-01', '1990-01-01',
       '1991-01-01', '1992-07-01', '1993-07-01', '1994-07-01',
       '1996-01-01', '1997-07-01', '1999-01-01', '2006-01-01',
       '2009-01-01', '2012-07-01', '2015-07-01'], 
      dtype='|S10')



In [17]:

    
(57204-mjdref)*86400









    Out[17]:





1435708800.0



In [18]:

    
(57754-mjdref)*86400









    Out[18]:





1483228800.0



In [19]:

    
# table_to_hdu doesn't seem to preserve the units of the column, so fix that
hdunew.columns.change_unit('LEAPSECS','s')
hdunew.columns









    Out[19]:





ColDefs(
    name = 'DATE'; format = '10A'
    name = 'TIME'; format = '15A'
    name = 'MJD'; format = 'D'
    name = 'SECONDS'; format = 'D'
    name = 'LEAPSECS'; format = 'D'; unit = 's'
)



In [29]:

    
def update_caldb_leapsec(NewLeapSecondDate, NewLeapSecond, updater="MFC", outdir=".", clobber=True):
    """
    Given the date of a new leapsecond (for example 2017-01-01T00:00:00)
    creates a new leapsecond file for transfer to the caldb

    NewLeapSecDate = date of new leap second in ISOT format YYYY-MM-DDTHH:MM:SS
    NewLeapSecond = amount of new leap second (usually 1.0)

    writes a new FITS file to the current working directory by default
    """
    from astropy.io import fits as pyfits
    from astropy.time import Time
    from astropy.table import Table
    import ftputil
    import time
    #
    # this block retrieves the lastest leapsecond file from /FTP/caldb/data/gen/bcf
    # based on the leapsecond file naming convention, "leapsec_mmddyy.fits"
    #
    LSdir = "FTP/caldb/data/gen/bcf/"
    host = ftputil.FTPHost('heasarc.gsfc.nasa.gov', "anonymous", "mcorcoran@usra.edu")
    genbcf = host.listdir(LSdir)  # get directory listing
    LeapsecFileList = [f for f in genbcf if 'leapsec' in f]
    LeapsecFileYear = [y.split("_")[1].split(".fits")[0][4:6] for y in LeapsecFileList]
    LeapsecFileYear = [('19' + y if int(y) > 50 else '20' + y) for y in LeapsecFileYear]
    LeapsecFileMonth = [m.split("_")[1].split(".fits")[0][2:4] for m in LeapsecFileList]
    LeapsecFileDay = [d.split("_")[1].split(".fits")[0][0:2] for d in LeapsecFileList]
    maxjd = 0.0
    for i in range(len(LeapsecFileList)):
        tiso = LeapsecFileYear[i] + "-" + LeapsecFileMonth[i] + "-" + LeapsecFileDay[i]
        fjd = Time(tiso).jd
        if fjd > maxjd:
            maxjd = fjd
            LatestLSF = LeapsecFileList[i]
    print("Latest Leapsecond File = {0}".format(LatestLSF))

    hdu = pyfits.open("http://heasarc.gsfc.nasa.gov/" + LSdir + "/" + LatestLSF) # open the file

    orig_header = hdu[1].header
    mjdref = orig_header['MJDREF']

    UpdateDate = time.strftime('%Y-%m-%d %H:%M:%S')
    #outfile = 'leapsec_' + (time.strftime("%d%m") + (time.strftime("%Y"))[2:]) + '.fits'
    outfile = 'leapsec_' + NewLeapSecondDate[8:10] + NewLeapSecondDate[5:7] + NewLeapSecondDate[2:4] + '.fits'


    #
    # create new row for new leapsecond information
    #
    newdate = NewLeapSecondDate.split('T')[0]
    newtime = NewLeapSecondDate.split('T')[1]
    newmjd = Time(NewLeapSecondDate, format='isot').mjd
    newsecs = (newmjd - mjdref) * 86400
    newLS = NewLeapSecond
    NewLeapSecondRow = [newdate, newtime, newmjd, newsecs, newLS]
    #
    # append new leapsecond to table and write to output file
    #
    tbdata = hdu[1].data
    t = Table(tbdata)  # convert hdu data to a python Table to add new row
    t.add_row(NewLeapSecondRow)  # add row of data
    hdunew = pyfits.table_to_hdu(t)  # convert table back to hdu (with minimal header)

    hdunew.columns.change_unit('SECONDS', 's') # table_to_hdu doesn't seem to preserve the Unit
    hdunew.columns.change_unit('LEAPSECS', 's') # table_to_hdu doesn't seem to preserve the Unit

    hdunew.header = orig_header  # use header from original file
    hdunew.header['COMMENT'] = UpdateDate+": "+updater+" ADDED "+NewLeapSecondDate+" LEAP SECOND"
    hdunew.header['HISTORY'] = "File modified by user "+updater+" on "+UpdateDate
    pyfits.writeto(outdir + "/" + outfile, hdunew.data, hdunew.header, clobber=clobber, checksum=True)
    return outfile

if __name__ == "__main__":
    file = update_caldb_leapsec('2017-01-01T00:00:00', 1.0)
    print file









    



Latest Leapsecond File = leapsec_010117.fits
Downloading http://heasarc.gsfc.nasa.gov/FTP/caldb/data/gen/bcf//leapsec_010117.fits [Done]
leapsec_010117.fits